#1. 数据分布
library(openxlsx)
###处理ADdata###
ADdata <- read.xlsx("E:\\class3\\13-刘睿佳-结果.xlsx")
rownames(ADdata)=ADdata[,1]   #将第一列作为行名
ADdata=ADdata[,-1]            #将第一列删除
a <- unlist(ADdata)         #把ADdata转化为向量 方法二：a <- as.matrix(ADdata)
c <- a[a != 0]        #去掉值为0的数据
c <- log(c)           #取对数
pdf("E:\\class5\\ADdatahist.pdf")  #将hist图存为PDF
hist(c,label=T,main="Histogram of ADdata",xlab="log(ADdata)",col=cm.colors(n))
dev.off()   #关闭PDF文件
###处理GSE67835数据###
GSEdata <- read.xlsx("E:\\class4\\13-刘睿佳-实践4结果.xlsx")
GSEdata <- GSEdata[,-c(1:4)]   #删除第1到4列
b <- unlist(GSEdata)        #把GSEdata转化为向量
d <- b[b != 0]        #去掉值为0的数据
d <- log(d)           #取对数
pdf("E:\\class5\\GSE67835hist.pdf")
hist(d,label=T,main="Histogram of GSE67835data",xlab="log(GSE67835)",col=cm.colors(n))
dev.off()


#2. 火山图
install.packages("ggplot2")
install.packages("ggrepel")
library("ggplot2")
library("ggrepel")
load("E:\\class5\\class5_volcano.RData")
data <- prostat
data$log10P = -log10(data[,3])

###数据处理###

data$Sig = ifelse(data$P < 0.05 &    #P<0.05以及FC>1.2的数据标记为'up',FC<1/1.2的数据标记为'down'，其余为'no'
                    abs(data$FC) > log2(1.2),  #abs绝对值
                  ifelse(data$FC <log2(1/1.2),'Down','Up'),'no')  
data = data.frame(data)
table(data$Sig) #查看数据统计情况

###绘制火山图###
p1 <- ggplot(data, aes(x =FC, y=log10P, colour=Sig)) + #x、y轴分别为log2FC和-log10P，颜色根据"Sig"
  geom_point(alpha=0.65, size=2) +  #点的透明度、大小
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757")) + xlim(c(-1.5, 1.5)) +  #调整点的颜色和x轴的取值范围
  geom_vline(xintercept=c(log2(1/1.2),log2(1.2)),lty=4,col="black",lwd=0.8) + #添加x轴辅助线,lty函数调整线的类型
  geom_hline(yintercept = -log10(0.05), lty=4,col="black",lwd=0.8) +  #添加y轴辅助线
  labs(x="log2 fold change", y="-log10 P-value") +  #x、y轴标签
  ggtitle("Volcano plot") + #标题
  theme_bw() + # 主题
  theme(plot.title = element_text(hjust = 0.5),
        legend.position="right", 
        legend.title = element_blank()
  ) 

###添加基因名标记###
p2 <- p1 + geom_text_repel(
  data = subset(data, data$P < 0.05 & abs(data$FC) >= log2(1.2)),
  aes(label = ID), size = 3,
  box.padding = unit(0.5, "lines"),
  point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE )

jpeg("E:\\class5\\火山图.jpg",width=1000,height=1000,quality = 100)  #将图输出jpg格式
p2
dev.off()

#3.对ADdata进行层次聚类分析,画树状图
###蛋白信息聚类分析###
ADdata.scaled <- scale(ADdata)     #数据标准化
d<-dist(ADdata.scaled)       #计算欧氏距离
windows()
fit1<-hclust(d,method = "average")
plot(fit1,hang = -1,cex=.8,main ="ADdata Protein Hclust")  #蛋白信息层次聚类树状图
pdf("E:\\class5\\ADdata蛋白信息层次聚类热图.pdf")
heatmap(as.matrix(d))  #用热图表示
dev.off()

###病人信息聚类分析###
ADdata.scaledt <- t(ADdata.scaled)
d2<-dist(ADdata.scaledt)       #计算欧氏距离
windows()
fit1<-hclust(d2,method = "average")
plot(fit1,hang = -1,cex=.8,main ="ADdata Patient Hclust")  #病人信息层次聚类树状图
pdf("E:\\class5\\ADdata病人信息层次聚类热图.pdf")
heatmap(as.matrix(d2))  #用热图表示
dev.off()
